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INTRODUCTION 


The  primary  objective  of  this  study  is  to  accurately  predict  the  amplitudes  and  wave 
forms  of  low  frequency  gravity  waves  in  the  atmosphere  and  the  ionosphere,  which  are 
produced  by  earthquakes  and  explosions.  Since  these  waves,  which  increase  in  amplitude 
with  altitude  because  of  the  decreasing  density  of  the  atmosphere  with  height,  will  produce 
relatively  large  fluctuations  in  the  electron  densities  in  the  ionosphere,  the  disturbance  can 
be  sensed  by  standard  electromagnetic  sounding  techniques.  Hence,  sensitive  monitoring  of 
seismically  produced  atmospheric  disturbances  can  be  accomplished  by  EM  sounding 
methods.  An  objective  of  this  study  is,  therefore,  to  produce  predictions  of  the  principal 
ionospheric  disturbances  to  be  expected  from  different  kinds  of  shallow  seismic  sources,  so 
that  their  characteristic  atmospheric  wave  signatures  can  be  used,  along  with  seismic 
methods,  to  identify  and  characterize  the  sources.  In  addition,  such  predictive  capabilities 
provide  a  framework  for  studying  the  properties  of  the  atmosphere  and  ionosphere  using 
gravity  wave  excitation  produced  by  large  seismic  sources. 

Near-surface  sources  can  produce  waves  in  the  atmosphere  having  two  principal 
components  at  long  distances  from  the  source,  an  acoustic  wave  with  relatively  high 
frequencies  and  a  gravity  wave  with  lower  frequencies.  Atmospheric  and  ionospheric 
disturbances  due  to  seismic  sources  have  been  observed  by  numerous  techniques  and 
experiments.  A  review  of  the  earlier  observations  is  given  by  BLANC  (1985),  where 
earthquakes  and  volcanic  eruptions  were  found  to  have  produced  atmospheric  and 
ionospheric  disturbances  both  near  to,  and  very  far  from,  their  epicenters  (DAVIES  and 
BAKER  1965).  In  several  early  studies,  nuclear  explosions  in  the  atmosphere  were  observed 
through  their  ionospheric  effects  (BENYON  and  JONES,  1962,  BAKER  and  DAVIES,  1968, 
FRANCIS,  1975).  There  have  also  been  recent  observations  of  acoustic-gravity  waves  in  the 
thermosphere  following  space  shuttle  ascents  (JACOBSON  and  CARLOS,  1994). 

Electromagnetic  sounding  of  the  ionosphere  after  surface  explosions  has  concentrated  on 
the  initial  acoustic  response  to  the  blast  wave  in  the  local  region  around  the  source  ( e.g . 
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BLANC  and  RICKEL,  1989).  However,  a  surface  chemical  explosion  of  5  kilotons  (ANFO) 
was  observed  to  excite  6  minute  oscillations  of  the  ionosphere  that  lasted  for  about  1  hour 
(JACOBSON,  CARLOS  and  BLANC,  1988).  KANAMORI  et  al„  (1994),  have  evaluated 
observations  for  several  major  volcanic  eruptions  and  have  found  5-minute  oscillations 
recorded  near  the  Mt.  St.  Helens  and  Krakatoa  volcanic  eruptions  with  amplitudes  of  0.3 
mbar  and  2  mbar  respectively.  Further,  larger  scale  ionospheric  motions,  produced  by  long 
period  gravity  waves  from  earthquakes  and  subsurface  nuclear  explosions,  are  likely  to  be 
responsible  for  the  observations  of  electromagnetic  signals  reported  by  GOKHBERG  et  al. , 
(1990). 

The  significance  of  the  5-minute  oscillations  is  that  this  period  closely  corresponds  to 
the  theoretical  Brunt-Vaisala  period;  that  is,  it  is  near  to  the  theoretical  minimum  period  for 
internal  gravity  wave  oscillations  in  the  atmosphere,  for  its  particular  density  decrease  with 
height.  Extensive  observations  of  the  distribution  of  gravity  waves,  both  in  frequency  and 
horizontal  wave  number  space,  indicate  that  the  observed  minimum  period  at  the  shorter 
wavelengths  is,  in  fact,  5  minutes  (MANSON,  1990).  Likewise,  the  observed  frequency 
distribution  of  gravity  waves  has  been  shown  to  be  a  consequence  of  both  wave  growth,  due 
to  decreasing  atmospheric  density  with  height,  and  wave  dissipation  (c.f.  EMAMURA  and 
OGAWA,  1995). 

Since  observations  of  zonal  winds  at  86  km  altitude  have  shown  that  the  background 
spectral  density  is  a  minimum  at  the  Brunt-Vaisala  frequency,  with  a  spectral  slope  of  about 
-5/3  for  frequencies  higher  than  those  in  the  tidal  band  (FRITTS,  1984),  it  follows  that  the  5 
to  10  minute  period  band  is  relatively  noiseless  and  so  gravity  waves  in  this  period  range 
should  be  more  readily  observable.  In  this  regard,  recent  observations  and  analysis  of 
nightglow  and  noctilucent  clouds  at  80  to  130  km  altitudes  has  allowed  the  direct 
observation  of  such  waves  (c.f.  GARDNER,  1995,  and  associated  papers).  In  addition,  since 
gravity  waves  refract  along  the  Earth's  surface,  they  are  observed  at  greater  horizontal 
distances  than  the  acoustic  wave  (FRANCIS,  1975).  Thus,  it  is  evident  from  the  standpoint 
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of  detectability  that  the  5  to  10  minute  period  gravity  wave  produced  by  a  seismic  source 
will  be  the  most  easily  measured  at  high  altitudes  near  the  source. 

As  noted  above,  it  has  been  demonstrated  that  wave  growth,  due  to  the  decrease  with 
altitude  of  the  density  of  the  atmosphere,  occurs  exponentially  with  a  lower  scale  height  of 
about  8  km,  and  it  is  also  known  that  this  growth  is  counteracted  by  damping  due  to  radiation 
and  molecular  diffusion  processes  below  the  ionosphere,  with  the  damping  magnitude  being 
dependent  on  the  frequency  and  wavelengths  of  the  internal  gravity  wave  ( IMAMURA  and 
OGAWA,  1995).  However,  gravity  waves  near  the  Brunt-Vaisala  frequency,  and  of  shorter 
wavelength,  are  less  affected  by  these  damping  processes  and  can  grow  to  a  saturation  point 
where  wave-breaking  with  viscous  damping  and  other  nonlinear  processes  can  occur. 
Observations  indicate  that  the  saturation  amplitudes  of  the  internal  gravity  waves,  occurring 
at  about  80  km,  are  approximately  10%  of  the  background  fluctuations  in  the  temperature 
and  pressure  fields  at  that  altitude,  with  the  gravity  wave  amplitude  decreasing  above  this 
altitude  (FRITTS,  1984).  Further,  the  dissipative  forces  are  found  to  be  such  that  short- 
duration  sources  produce  far-field  gravity  waves  with  wavetrains  from  1  to  2  hours  in 
duration.  In  the  non-linear  modeling  developed  in  the  present  study,  we  will  incorporate 
specific  damping  mechanisms  and  compare  predicted  results  with  observations  in  order  to 
confirm  that  the  modeling  reproduces  observations. 

Based  on  partial  confirmation  of  the  non-linear  numerical  modeling  by  these  limited 
observational  results  and  verified,  albeit  approximate,  analytical  results,  we  will  use 
predictions  of  amplitudes  and  wave  characteristics  in  the  time  domain  to  evaluate  the 
possibilities  for  observational  studies  using  particular  atmospheric  and  ionospheric  gravity 
wave  detection  methods.  In  this  regard,  we  compare  correlations  between  the  non-linear 
gravity  wave  predictions  of  this  study  and  the  GPS  detections  of  ionospheric  perturbations 
obtained  by  Calais  and  Minster  (1995)  following  the  large  January  17,  1994  Northridge 
earthquake  in  order  to  try  to  establish  a  framework  for  further  investigations. 
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CONSERVATION  RELATIONS  FOR  ATMOSPHERIC  MOTIONS 


We  use  the  primitive  system  where  the  continuity  equations  are  based  on  the  values  of 
the  continuum  fields  at  particular  points  in  space  and  time.  Conservation  of  mass  is  expressed 
as  : 


f  +  ^(pvP  =  0  (1) 

where  p  is  density  and  Vj  is  velocity  in  the  Xj  direction.  Conservation  of  momentum  is 
expressed  as  : 


dvj  dv; 

+  V^J 


=  pXi  +  -£Pu  +  Fi(d) 


(2) 


where  X;  are  external  forces  such  as  gravity,  and  Py  is  the  generalized  stress,  with: 


P.j  =  “  P-5y  +  2n-ey  -  2/3p-6ij  eki 


(3) 


Here.  e;j  =  1/2 


is  the  strain  rate  and  p  the  viscosity.  An  effective  force  term  Fj(d) 


[dxj  3x,j 

is  included  to  account  for  microscale  momentum  transfer  that  damps  the  mean  velocity  field  v; 
appearing  in  (2).  This  force  will  be  referred  to  as  a  “drag  force"  and  is  usually  taken  to  be 
proportional  to  velocity,  with  various  phenomena  such  as  turbulence  and/or  ion  drag  contribut¬ 
ing  in  different  altitude  regions. 

Conservation  of  energy  gives  the  usual  continuity  equation  governing  temperature: 

p i  ^ + -4  <«■  *>  -  $  -  P'S + *M> * 

where  T  is  temperature.  Cy  is  the  specific  heat  at  constant  volume,  K  is  the  thermal  conduc¬ 
tivity  and  <l>(d)  the  dissipation  function  involving  drag  and  <t>(w)  the  dissipation  function  involv¬ 
ing  viscous  loss  terms. 

The  equation  of  state  for  the  atmospheric  gas  is  taken  to  be  ideal  i.e. 


_  _ 
P  =  —  p-T 
m 


(5) 
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where  kg  is  Boltzmann's  Constant  and  m  is  the  mean  molecular  weight. 

In  the  fundamental  equations,  (1)  through  (5),  the  dependent  and  independent  variables 
are  now  normalized  with  respect  to  typical  values.  For  an  ambient  atmosphere  with  exponen¬ 
tial  decay  of  density  with  height,  distances  arc  normalized  through  the  scale  height,  H,  which, 
at  the  surface,  is  approximately  8400  metres.  Velocities  are  normalized  with  respect  to  c,,  the 
sound  velocity  of  air  at  the  earth’s  surface.  Thus,  the  Mach  Number  is  just  v/c,  where  v  is  a 
typical  velocity  of  the  flow.  Similarly  density,  pressure  and  temperature  are  normalized  to 
ambient  surface  values,  and  the  independent  time  variable,  t,  is  normalized  by  (H/cJ.  This  pro¬ 
duces  a  set  of  dimensionless  ratios  related  to  the  Mach  and  other  well-known  dimensionless 
numbers. 


Thus,  for  the  continuity  of  mass  we  get.  as  before. 


dp  9  ,  . 

a.  =  <« 

where  the  new  variables  are  now  normalized  and  given  the  same  symbol  as  the  original  vari¬ 
ables.  Incorporating  gravity  as  an  external  force,  the  momentum  conservation  equation,  (2). 
becomes  upon  normalization  : 


9v;  dVj  a2  dP  A4  A, 

’  =  ”  Vj'  97  -G*  S(z)'e*  “  ~ +  —•%  +  —  Fi 


at 


P  dxi 


(7) 


where  G,  =  (g^H/c,2)  is  a  combination  of  Mach  Number  and  Froude  Number  and  the  dimen¬ 
sionless  number  A2  =  p^/Cp^-c,2)  combines  the  Euler  and  Mach  Numbers.  A4  =  p^/PjC,  meas¬ 
ures  the  ratio  of  viscous  to  inertial  forces  and  combines  the  usual  Reynolds  and  Mach  Numbers 
and  'F;  is  the  normalized  viscous  force,  given  by  the  spatial  gradient  of  the  viscous  terms  in  the 
generalised  stress  equation  (3)  and  where  |i  is  normalized  to  |ts  ,  the  viscosity  at  the  surface. 
In  the  atmosphere,  p  is  usually  taken  to  be  constant  for  the  molecular  viscosity.  The  parameter 
A,  =  H/(pscs2)  normalises  the  drag  force  whose  particular  dependence  will  be  discussed  later. 
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In  the  case  of  the  conservation  of  energy,  the  normalized  equation  is 


5T 

3t 


+  Ag- 


a^r 

ax/ 


-A«- 


a^ 

ax: 


+  a5  <i>(v)  +  a7  <j>(d) 


(8) 


where 


As  =  H, 


kB  C»  _  a  A< 
m-Cy  (p,H)  ^  A2 


Ats  - 


m-Cy 

A7  Cj/fCyPjHX,) 

T. 


Ag  =  A*  •  K 


(P.C.H) 

Here  the  thermal  conductivity  K  is,  as  usual,  taken  to  be  constant  for  the  atmosphere.  The 
equation  of  state,  on  normalization,  is 


p  =  p*T/m(z) 


(9) 


since  the  ideal  equation  of  state  at  the  surface  is  p,  =  —  •  p,  •  T$  where  m,  is  the  surface 
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value  of  mean  molecular  weight  and  m(z)  is  the  height-dependent  normalized  value. 


Because  an  acoustic-gravity  wave  is  a  perturbation  on  the  existing  atmospheric  structure, 
we  shall  decompose  the  thermodynamic  variables  into  two  components;  the  background  sta¬ 
tionary  component  and  the  perturbation.  However,  we  do  not  assume  that  the  perturbation  is 
small  in  either  temporal  or  spatial  distribution  as  the  gravity  wave  can  grow  with  altitude. 
Thus,  expressing  the  normalised  density  as  p  =  p0  +  Pi.  where  0  indicates  the  stationary  state 
and  1  indicates  the  perturbation,  we  get  for  the  conservation  of  mass  : 

ir  “  -  i;[<p»  * Pi,v'1  <,0) 

Similarly,  for  the  velocity  equation  (7).  we  get : 


at 


=  —  V;- 


avj  __Gfp, 

9X: 


a2  ap  i 

•g(zK - — -T-  + 


•4':  + 


;(«1) 


-j  Po  +  Pi  *  Po  +  Pi  ^xi  Po  +  Pi  Po  +  Pi 
where  P[  is  the  perturbed  pressure  component  The  temperature  equation  (8)  becomes  : 


(ID 
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(12) 


3T, 


9(T0  +  T,)  +  Aj.  JL<x0  +  Ti)  _  +  Aj.<j>cv)  +  A?.4>W) 


9xi 


dXi 


*j  UAj  — j 

where  T,  is  the  perturbed  temperature  component.  The  equation  of  state  (9)  becomes  : 


Pi  =  (PoT,  +  PiT0  +  p,T,)/m(z) 


(13) 


The  introduction  of  internal  sources  and  sink's  of  energy  and  momentum  is  intended  to 
account  for  processes  which  are  poorly  understood  theoretically  and  are  represented  by  parame¬ 
terization  (VOLLAND,  1988).  In  the  atmosphere  below  the  thermopause,  the  only  dissipative 
mechanisms  arc  molecular  viscosity,  thermal  conductivity  and  interaction  with  the  innate  tur¬ 
bulence.  The  kinematic  viscosity  increases  like  the  density  decreases  with  height  in  the  atmo¬ 
sphere.  However,  neither  it  nor  the  molecular  thermal  conductivity  immediately  affect  the 
momentum  and  energy  of  these  gravity  waves  below  the  thermosphere,  as  it  takes  a  time  of  the 
order  of  days  to  damp  these  large  wavelength  systems  (PITTEWAY  and  HINES,  1963).  Thus, 
we  ignore  the  molecular  viscosity  term  %  in  equation  (1 1)  and  the  <J>(V)  term  in  equation  (12). 
Above  the  thermopause,  it  is  expected  that  dissipation  due  to  ion  drag  becomes  important 
(MIESEN  et  at..  1989). 

Turbulence,  at  altitudes  below  the  thermopause,  is  considered  to  be  the  most  important 
drag  force  (DEARDORFF,  1985).  In  order  to  account  more  accurately  for  the  inherent  sub¬ 
grid  scale  turbulence  in  the  atmosphere  below  the  thermopause,  the  flow  variables  at  a  point 
can  be  decomposed  into  a  mean  flow  and  a  perturbed  turbulent  flow  (DEARDORFF.  1985). 
In  the  momentum  continuity  equation,  additional  components  are  thus  obtained  for  the  general¬ 
ized  stress,  which  are  mainly  interaction  terms  between  the  mean  and  perturbed  densities  and 
velocities,  termed  the  Reynold’s  stresses,  which  represent  the  interaction  of  the  mean  flow  with 
the  sub-grid  scale  turbulence.  These  additional  stresses  have  been  approximated  in  various 
phenomenological  approaches,  including  the  introduction  by  Boussinesq  of  the  concept  of  eddy 
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viscosities  in  order  to  use  the  Newtonian  equations  with  the  usual,  but  much  larger,  viscosity 
term.  However,  this  equivalent  viscosity  approach  has  been  shown  to  be  unable  to 
effectively  model  dissipation  due  to  turbulence,  and  so  other  methods  using  Rayleigh  friction 
have  been  invoked  (VOLLAND,  1988). 

We  will,  therefore,  model  the  effects  of  turbulence  on  the  gravity  wave  momentum 

through  use  of  drag  forces  that  are  proportional  to  velocity.  We  incorporate  these  diverse 
sources  of  drag  and  dissipation  into  the  force  term  Fjd)  that  acts  to  reduce  the  mean 

velocity.  That  is,  in  eq:  (7)  we  use: 

F-d)  =  -k(d)  •  Vj  (14) 

where  kW)  is  a  spatially  dependent  drag  coefficient.  Different  drag  components,  due  to 
differing  mechanisms,  will  be  operable  at  various  altitudes  and  under  various  wind  and  ion 
conditions,  so  that  k  (d)  can  be  explicitly  altitude  dependent.  The  effect  of  the  dissipative 
mechanisms  on  the  temperature  equations  are  neglected,  as  we  invoke  the  argument  that 
energy  loss  from  the  large-scale  gravity  waves  through  interactions  with  large-scale  turbulent 
motions  takes  a  time  of  the  order  of  days  to  decay  into  molecular  size  eddies  and  thereby 
input  thermal  energy  into  the  system  (PITTEWAY  and  HINES,  1963).  Thus,  the  <b(d)  term 
in  eq:  (12)  will  be  neglected  compared  to  other  terms,  since  short  term  predictions  are  of 
interest. 

The  boundary  conditions  that  apply  express  conservation  of  mass  and  momentum  across 
surfaces  of  material  discontinuity.  They  are  simply  expressed  as  the  continuity  of  the  normal 
component  of  particle  velocity,  for  conservation  of  mass,  and  the  continuity  of  tractions,  for 


the  conservation  of  momentum.  That  is: 


where  r  denotes  the  position  vector,  n  is  the  normal  to  the  surface  S,  and  the  double  bracket 
notation  in  (15)  is  used  to  represent  the  difference  in  the  bracketed  quantity  across  the 
surface  S. 

A  seismic  source,  producing  waves  in  the  solid  medium  below  the  planetary 
atmosphere-lithosphere  boundary  (S0),  produces  atmospheric  excitation  through  the 
continuity  equations  in  (15).  From  these  conditions  we  see  that  the  vertical  component  of 
the  solid  medium  particle  velocity  is  imparted  to  the  atmospheric  medium,  while  the  other 
coupling  is  that  of  the  atmospheric  pressure  to  the  solid  medium  traction  on  S0-  To  first 
order  the  pressure  changes  in  the  atmosphere  due  to  the  traction  variations  produced  by 
seismic  waves  are  neglected  and  the  boundary  is  treated  as  traction  free.  In  this  case,  only 
the  velocity  condition  in  (15)  produces  time  dependent  momentum  changes  at  the  lower 
boundary  of  the  atmosphere. 

Thus,  the  atmospheric  coupling  due  to  a  seismic  source  in  the  earth  is  quite  predictable 
when  the  seismic  wave  field  in  the  solid  can  be  specified  either  numerically  or  analytically. 
Representations  of  wave  fields  from  seismic  sources  of  various  kinds  (e.g.  earthquakes, 
chemical  and  nuclear  explosions)  can  be  made  with  accuracy,  so  that  predictions  of 
atmospheric  excitation  using  the  velocity  boundary  condition  in  (15)  is  quire  straightforward. 
Here,  the  direct  wave  field  from  the  seismic  source  plays  the  role  of  a  forcing  function  on  the 
atmosphere  along  the  boundary  at  S0,  and  will  give  rise  to  low  frequency,  long  wavelength, 
atmospheric  motions,  including  net  flow  components,  as  well  as  high  frequency  acoustic 
waves.  In  the  case  of  earthquakes,  and  underground  nuclear  tests,  which  release  most  of 
their  energy  at  depths  well  below  S0,  the  excitation  of  high  frequency  acoustic  waves  can  be 
expected  to  be  negligible  compared  to  transport  effects  and  the  boundary  movement 
excitation  of  the  much  lower  frequency  gravity  waves. 

The  atmospheric  motions  are  also  the  cause  of  ionospheric  charged  particle  motions 
which  can  be  approximated  as  following  the  motions  of  the  neutral  gas  (KELLEY  and 
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HYSELL,  1991).  The  basic  conservation  law  for  charged  particles,  assuming  no  creation  or 


destruction,  is: 


3Na  a(Na-v“) 

dt  dx. 


where  Na  is  the  number  of  particles  of  type  a  and  v“  is  their  velocity  in  the  j'th  direction. 

The  particles  of  interest  here  are  electrons,  and  the  initial  concentration  of  this  charged 
particle  is  taken  to  be  time-independent  with  only  a  vertical  functional  dependence,  N(z). 
Since  we  are  concerned  with  one  particle  type  the  subscript  a  will  be  dropped  in  the 
following.  Assuming  only  small  changes  in  this  concentration,  the  time  dependence  can  be 
found  from  integrating  equation  (16)  over  time.  In  this  approximation  the  concentration 


change  becomes: 


5N(z,t)  =  -^^-j  vz(t')  dt'-N(z)- j  ^pdt'  (17) 

to  t0  J 

The  first  term  in  (17)  is  the  concentration  change  due  to  the  displacement  of  the  ionospheric 
layer,  while  the  second  term  arises  as  a  result  of  compression  or  rarefaction  and  is  the 
dominant  term  when  dealing  with  processes  involving  characteristic  dimensions  smaller  than 
the  width  of  the  layer.  The  velocity  of  the  charged  particles  is  usually  taken  to  be  of  the 
order  of  that  of  the  neutral  gas  and  this  approximation  will  be  used  in  the  calculations  for 
electron  density  charges. 

Therefore  the  neutral  gas  velocity,  which  is  computed  using  equ:  (10)  -  (13)  is  used  in 
eq:  (17)  to  predict  the  electron  density  changes  in  the  ionosphere.  In  this  latter  computation, 
the  initial  concentration  of  electrons  is  taken  to  be  that  of  a  Chapman  distribution.  This 
distribution  has  a  maximum  density  at  345  km  height  and  decreases  rapidly  below  about  90 
km  ,with  the  functional  dependence  on  height  defined  by: 


N(z)  =  Nc  -  exp 


where  ^=(z-hc)/H,  hc  is  345  km,  H  is  65  km  and  Nc  is  a  normalizing  value  equal  to  N  at  the 
altitude  hc 
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Linear  Approximations  for  Gravity  and  Internal  Waves 

Predictions  of  important  features  of  the  low  frequency  gravity  and  internal  waves 
generated  by  surface  atmospheric  sources  have  been  made  using  linearized  approximations 
to  the  eq:  in  (10)  through  (13).  Many  of  these  predictions  have  been  observationally 
verified,  at  least  in  the  approximate  sense  of  the  theory  itself.  Therefore,  the  verified 
analytical  results  from  linear  approximations  are  important  for  physical  insights  and  as 
criteria  in  evaluating  results  from  non-linear  numerical  computations. 

Although  the  full  set  of  flow  equations  do  not  have  analytical  solutions,  it  has  been 
shown  that  even  weak  non-linear  analyses  can  yield  solitary  waves  (MDESEN,  1992)  and 
complex  interactions  ( KLOSTERMYER,  1978).  Linearization  of  the  fundamental 
conservation  equations  is  usually  performed  by  removing  terms  that  are  of  second  order  in 
particular  physical  situations.  This  approximation  reduces  the  equ:  (10)  -  (13)  to: 


9xj 

3vi  ,  .  3P, 


+  Y,+F 


'(d) 


(19) 

(20) 


which  apply  only  for  flows  and  waves  with  quite  low  particle  velocities  and  velocity 
gradients. 

Internal  gravity  waves  in  density-stratified  fluids  are  completely  different  from  the 
familiar  acoustic  waves,  as  they  are  anisotropic  and  dispersive.  After  linearization,  they  are 
governed  by  hyperbolic  differential  equations,  rather  than  the  elliptical  equations  obtained  in 
the  linearized  acoustic  case.  Thus,  the  gravity  waves  are  anisotropic  and  dispersive  and  do 
not  obey  Fermat's  principle  (BARCILON  and  BLEISTEIN,  1969).  Fourier  decomposition  of 
the  linearized  perturbation  equations,  with  phase  and  group  velocities  assuming  paramount 
importance,  has  been  successful  in  analyzing  these  internal  waves  (LIGHTHILL,  1978). 
Here,  the  emphasis  has  been  on  initial  perturbation  problems  and  oscillating,  or  uniformly 
moving,  sources.  In  the  simplest  case,  the  fluid  is  assumed  unbounded,  dissipative  forces  are 
ignored,  the  buoyancy  (Brunt-Vaisala)  frequency  is  assumed  constant  and  the  internal  waves 
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are  three-dimensional.  Further,  the  Boussinesq  approximation  is  assumed  and  this  implies 
the  neglect  of  the  inertial  effects  of  density  variations  compared  with  the  buoyancy  forces 
they  create.  Lighthill  (1978)  also  argued,  however,  that  when  non-Bousinesq  effects  are 
taken  into  account,  compressibility  is  important  and  the  internal  waves  become  acoustic - 
gravity  waves. 

In  a  stratified  fluid,  where  the  undisturbed  density  p0  varies  exponentially  with  height  z 
according  to 

Po(z)  =  p00e_pz 

the  buoyancy  frequency  N  is  given,  in  the  linear  approximation,  by 

N  =  (g(3)1/2 

For  the  simplest  case,  the  dispersion  relation  for  plane  monochromatic,  Boussineq  type 
internal  waves  (LIGHTHILL,  1978)  is  given  by 

(0  =  Nkh  /  k 

where  k  is  the  modulus  of  the  wavevector  and  kh  is  the  modulus  of  its  horizontal  projection. 
This  implies  waves  with  a  frequency  co<N,  for  an  arbitrary  wavelength  X  ,and  an  inclination 
of  the  planes  of  constant  phase  to  the  vertical  of  0  =  arccos  (CD  /  N).  The  phase  velocity 
with  which  these  planes  move  is  perpendicular  to  the  group  velocity  and,  consequently,  the 
fluid  particles  move  along  straight-line  paths  parallel  to  the  wavecrests.  In  the  non- 
Boussinesq  case,  the  dispersion  relation  becomes  (VOISIN,  1991): 

CD  =  Nkh  /  (k2  +  0.25p2)1/2 

Thus,  the  wavenumber  surface  is  no  longer  a  cone  but  a  hyperboloid  surface  of  revolution. 
The  group  velocity  points  along  this  surface  normal,  but  it  is  no  longer  perpendicular  to  the 
wavevector,  and  the  trajectories  of  fluid  particles  now  become  ellipses.  A  point 
monochromatic  source  thus  radiates  non-Boussinesq  internal  waves  into  the  total  region 
where  |cos0|<  CD  /  N,  with  a  group  velocity  that  is  a  maximum  at  a  particular  angle  in  this 

cone. 
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The  more  recent  linear  approach  by  VOISIN  (1991)  is  based  on  a  Green’s  function 
method,  where  the  internal  wave  field  is  studied  using  a  number  of  different  Green's 
functions  applicable  under  various  specific  conditions.  Voisin  shows  that,  in  accordance 
with  the  group  velocity  approach  used  by  Lighthill,  internal  waves  propagate  for  frequencies 
in  the  range  N|cos0|<|(o|<N,  while  outside  this  band  they  are  evanescent  with  an 
exponential  decay  with  distance  from  the  source.  Further,  Voisin  determines  that  impulsive 
internal  waves  are  made  up  of  the  combination  of  gravity  waves  and  buoyancy  oscillations  of 
the  fluid;  with  the  latter  described  by  DICKINSON  (1969). 

In  order  to  understand  the  separation  of  internal  waves  into  gravity  waves  and  buoyancy 
oscillations,  asymptotic  approaches  are  used.  For  a  point  source,  at  small  times,  the 
Boussinesq  fluid  motion  essentially  ignores  its  internal  stratification  and  its  motion  is 
irrotational  (BATCHELOR,  1967).  At  large  times  such  that  Nt  »1,  gravity  waves  and 
buoyancy  oscillations  become  separated,  with  gravity  waves  having  the  phase  expected  from 
group  velocity  analysis  while  buoyancy  oscillations  are  present  but  do  not  propagate.  Thus, 
at  large  times,  non-Boussinesq  effects  are  such  that  the  internal  waves  gradually  build  up 
(TOLSTOY,  1973)  and  eventually  these  internal  waves  separate  into  two  components  which, 
in  the  limit  of  large  time,  are  identified  as  gravity  waves  and  buoyancy  oscillations. 
Buoyancy  oscillations  in  this  case  are  waves  which,  unlike  gravity  waves,  consist  of  both 
propagating  and  evanescent  internal  waves,  the  later  decaying  exponentially  in  time.  For 
large  times,  the  gravity  and  buoyancy  waves  split  and  increasingly  separate  and  eventually 
lose  their  non-Boussinesq  character.  In  this  limit,  the  Boussinesq  gravity  waves  are 
(approximately)  plane  propagating  internal  waves  of  frequency  N|cos0|  and  wavelength 
2ra7  NtsinO.  On  the  other  hand,  the  Boussinesq  buoyancy  waves  are  radial  oscillations  of 
frequency  N  and  are  present  everywhere  in  the  fluid  with  a  wavevector  that  is  horizontal.  In 
the  non-Boussinesq  far-field  ((3r»l),  VOISIN  (1991)  has  shown  that  a  vertical  point  force 
(F0)  at  the  ground  level  generates  a  pressure  field  which  decays  as  r  "2  and  is  proportional  to 
F  o- 
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Numerical  Modeling  of  Non-Linear  Atmospheric  and  Ionospheric  Flows  and  Waves 

The  set  of  non-linear  partial  differential  equations  governing  the  atmospheric  flows  are 
converted  to  a  corresponding  set  of  finite  difference  equations  in  order  to  perform  explicit 
computer  integration  in  time  and  space.  By  examining  exact  solutions  to  linear  and  non¬ 
linear  partial  differential  equations  and  their  corresponding  difference  equations,  two  rules  of 
modeling  have  been  formulated  (MICKENS,  1988).  In  order  to  prevent  "ghost"  solutions 
(numerical  instabilities),  the  first  rule  is  that  the  order  of  the  finite-difference  scheme  should 
be  equal  to  the  order  of  the  differential  equation.  The  second  rule  applies  to  non-linear 
components  of  the  differential  equation  and  requires  such  non-linear  terms  to  be  treated  non- 
locally  on  the  lattice  for  stability  (MICKENS,  1988). 

For  example,  the  difference  equation  used  for  the  velocity  equation  in  the  x  direction 
can  be  formulated  as: 


chi  du  du  du  _ 
— +  u—  =  -v-  —  -w—  +  F 
dt  ox  dy  dz 


(21) 


where  F  represents  the  remaining  terms  on  the  right-hand  side  of  the  x  component  of  eq: 
(11).  Upwind  differencing  is  used  for  first  order  spatial  gradients  in  the  corresponding 
difference  equation.  However,  where  the  velocity  operates  on  its  own  velocity  gradient,  such 
non-linear  terms  are  treated  non-locally  in  a  fashion  similar  to  that  which  is  necessary  when 
integrating  the  Korteweg-deVries  equation  for  multiple  soliton  solutions  (CRANDALL, 
1991).  The  difference  equation  for  (21)  is  thus: 


Uijk  Uijk  n+1  Uijk  Ui-ljk  _  n  Uijk  Uij-lk  n  Uijk  Uijk-1  pD 

. . XT  u,)l  _ ta  v«‘  Ty  «l  S  Fi|k 


(22) 


Hence,  the  new  value  of  the  velocity  component,  u,  multiplies  the  x  component  of  the 
spatial  gradient  composed  of  old  u  values,  and  so  results  in  a  non-linear  integration  algorithm 
given  by: 


-  v5k Kk  -  Uy_u) 7^  -  w,5k(uSk  —  U^)-^ +  AtF,1 

n+1  _  _ Ay _ AZ 


n 

ijk 


Uijk  = 


l  +  (u«k  -u“_ijk)At/  Ax 


(23) 
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Because  we  are  computing  density  changes  rather  than  the  density  itself,  which  is 
always  positive,  it  is  not  necessary  to  use  an  algorithm  which  keeps  the  density  from  going 
negative,  a  restriction  which  can  cause  much  difficulty  (ORAN  and  BORIS,  1987).  The 
updated  variable  is  projected  not  from  just  the  old  dependent  variable,  a  process  that  is 
inherently  unstable,  but  from  a  distributed  smoothed  average  of  the  variable  at  locations 
surrounding  the  specific  spatial  location.  Such  a  smoothing  method  brings  stability  to  the 
differencing  scheme.  However,  the  attendant  numerical  diffusion  is  minimized  through  anti- 
diffusion  techniques  (ORAN  and  BORIS,  1987)  incorporated  into  the  integration  algorithm. 

The  second  order  derivatives  in  the  viscosity  and  thermal  conductivity  terms  are  also 
modeled  by  finite  differences  taken  at  the  surrounding  spatial  locations.  In  the  explicit 
integration  scheme,  the  updated  flow  velocities,  temperature  and  density  perturbations  are 
obtained  via  their  continuity  equations,  while  the  updated  pressure  perturbation  is  obtained 
from  insertion  of  the  updated  density  and  temperature  into  the  ideal  gas  equation. 

Aside  from  the  various  differential  equations,  there  are  at  least  three  boundary 
conditions  important  to  the  modeling  of  fluid  flows.  At  the  bottom  boundary  a  vertical 
velocity  variation  along  the  boundary  is  prescribed  by  the  seismic  wave  field,  as  previously 
discussed.  (Because  of  the  presence  of  a  lower  boundary  layer  above  a  complex  topography, 
horizontal  velocities  are  not  constrained  to  be  zero  but,  instead,  constant  velocity  and 
density  gradients  are  assumed  in  the  vertical  direction.)  On  the  other  hand,  the  topmost 
boundary  should  mimic  the  conditions  for  an  open  atmosphere,  with  specific  considerations 
for  buoyancy  and  field  gradients.  We  have  examined  various  options  including  fixing 
velocities,  densities  and  their  gradients.  In  the  end,  after  many  tests,  we  have  adopted  a 
general  open  flow  boundary  condition,  which  we  also  use  for  the  artificial  side  boundaries, 
where  all  normal  gradients  are  constant  at  these  boundaries.  However,  if  applied  to  the 
temperature  variable,  this  would  not  permit  heat  flow  through  the  open  boundary.  Therefore, 
for  temperature,  the  second-order  normal  derivative  is  made  constant. 
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Figure  1.  Cylindrically  symmetric  two-dimensional  solution  for  the  pressure  field  due  to  a 
constant  velocity  air  source  at  the  base  of  the  atmosphere. 
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In  general,  a  sufficiently  large  source  near  the  surface  will  cause  vertical  velocity  and 
density  variations  in  the  atmosphere  which  move  upward  and  are  initially  observed  as  vortex 
ring  thermals.  Here,  the  full  non-linear  equations  are  needed,  since  the  advective  terms 
become  of  importance  inside  this  vortex  ring  (TURNER,  1973).  Experiments  (SCORER, 
1957)  show  that  a  vortex-like  circulation  is  superimposed  on  the  initial  general  vertical 
motion  during  the  expansion  period  of  the  thermal,  whose  shape  is  slightly  oblate  spheroidal. 
Mixing  takes  place  with  the  fluid  ahead  of  the  advancing  front  of  the  thermal  and,  as  a 
result,  the  volume  of  the  thermal  continuously  grows  but  the  velocity  eventually  decreases. 

Numerical  simulations  predicting  the  motion  and  structure  of  these  thermals  can  be  used 
as  a  test  of  modeling  accuracy.  In  this  regard,  a  vertical  cross-section  of  a  typical  example 
is  shown  in  Figure  1.  As  a  function  of  time  the  transients  propagate  upward,  with  advective 
motion  occurring  both  vertically  and  horizontally.  A  sequence  of  spiral  circulation  patterns 
develops  through  asymmetric  screw  modes  across  the  horizontal  cross-section,  as  previously 
argued  by  NYGREN  et  al.  (1984),  and  energy  and  momentum  are  circulated  through  the 
plume  by  traveling  waves  through  the  circulation  pattern,  as  seen  by  WEIL  (1988).  Similar 
effects  have  been  observed  in  the  real  atmosphere  when  thermals  and  plumes  propagate 
upward  with  similar  asymmetric  shearing  flows  (KUETTNER  et  al.  1987).  At  the  neutral 
altitude,  the  buoyancy  is  zero  in  a  density  decreasing  atmosphere,  and  these  thermals  break 
down  into  internal  gravity  waves  through  oscillations  around  the  neutral  elevation. 
Consequently,  the  non-linear  modeling  approach  used  here  produces  results  that  agree  with 
observations  and  with  results  of  other  lines  of  analysis  as  well. 

Seismic  Source  Effects 

We  are  concerned  with  a  variety  of  near-surface  seismic  sources  that  produce  varying 
degrees  of  movement  of  the  earth-atmosphere  boundary  which,  in  turn,  result  in  the 
excitation  of  atmospheric  disturbances.  In  particular,  earthquakes  can  produce  large,  near- 
field,  displacements  of  the  surface  of  the  earth  surrounding  the  epicenter  of  the  event. 
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Further,  shallow  earthquakes  excite  high-amplitude  surface  waves  which  propagate  to  large 
distances.  Because  of  the  large  surface  displacements  associated  with  these  waves  and  the 
near-source  boundary  movements,  low-frequency  disturbances  in  the  atmosphere  are 
produced  which  are  quite  large. 

Underground  nuclear  explosions  also  produce  a  piston-like  motion  of  the  earth’s  surface 
directly  above  the  source  as  well  as  large  seismic  surface  waves  at  larger  distances. 
Specifically,  tamped  nuclear  explosions,  with  the  device  in  near  proximity  to  the  surrounding 
rock,  produce  high  particle  velocities  at  the  surface  boundary  and  strong  atmospheric 
coupling  while  "decoupled"  explosions,  which  are  detonated  in  large  cavities,  produce  much 
lower  boundary  velocities  and  little  atmospheric  effects.  Chemical  explosions,  for  industrial 
purposes,  produce  large  surface  motions  even  when  the  yield  is  small  since  they  are  so  near 
the  surface. 

Therefore,  all  seismic  sources  produce  surface  motions,  with  the  vertical  velocity 
component  of  the  boundary  surface  acting  as  a  piston  on  the  local  atmosphere.  In  this  case 
compressive  motions  decrease  the  buoyancy  of  the  local  atmosphere  while  downward  surface 
motions  expand  the  local  atmosphere  and  thereby  produce  increased  buoyancy.  In  order  to 
model  these  sources  in  numerical  codes,  it  is  necessary  to  model  the  time  dependent  particle 
velocity  and  stress  at  the  free  surface  of  the  earth  in  the  vicinity  of  the  source  epicenter.  As 
indicated  in  eq:  (15),  the  boundary  conditions  require  continuity  of  the  velocity  component 
normal  to  the  surface  and  continuity  of  the  tractions  at  the  surface.  Consistent  with  the  usual 
seismic  theory  approximations,  the  shear  tractions  in  the  atmosphere  are  taken  to  be 
negligible.  Likewise,  the  pressure  fluctuations  in  the  atmosphere  due  to  dynamic  traction 
variations  in  the  solid  are  taken  to  be  of  second  order  relative  to  the  pressure  variations 
caused  by  the  seismically  driven  vertical  movements  of  the  air-solid  boundary.  Thus,  the 
traction  continuity  requirement  is  automatically  satisfied,  to  first  order  at  least,  and  the 
remaining  condition  that  produces  first  order  coupling  effects  in  the  atmosphere  is  that  given 
by  the  continuity  of  the  normal  component  of  particle  velocity  at  the  boundary. 
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Earthquakes  produce  large  concentrated  movements  along  linear  faults  with  lengths 
which  typically  extend  from  tens  to  hundreds  of  kilometers  for  moderate  to  very  large  events. 
Near  epicenters  and  near  the  fault  trace  of  shallow  earthquakes,  ground  surface  accelerations 
significantly  greater  than  the  acceleration  of  gravity  at  the  surface  are  commonly  observed. 
Thus,  there  is  the  potential  for  large  boundary  velocities  and  strong  atmospheric  coupling. 

In  this  regard,  strike-slip  faulting  will  produce  no  vertical  motions  at  a  flat  free  surface, 
but  when  appreciable  topography  exists  such  motions  can  initiate  upward  flows  and  waves  in 
the  atmosphere.  On  the  other  hand,  shallow  dip-slip  and  thrust  faulting  always  will  produce 
large  vertical  motions  of  the  earth's  surface,  with  opposite  sides  of  the  fault  moving  in 
opposite  directions.  Thus,  earthquakes,  even  at  depths  that  produce  no  surface  breaks,  can 
produce  appreciable  near-field  surface  velocities  and  significant  static  surface  displacements 
close  to  the  epicenter,  as  well  as  large  low-frequency  surface  waves  propagating  thousands 
of  kilometers  from  the  source.  Both  near-source  boundary  movements  and  the  far  field 
surface  waves  have  been  observed  to  excite  atmospheric  waves  (DAVIES  and  BAKER, 
1965). 

Large  mining  explosions  also  produce  boundary  movements  with  an  areal-temporal 
distribution;  in  this  case  caused  by  ripple  firing  and  the  finite  time  delay  between  the 
sequence  of  separate  explosions  that  are  commonly  used.  These  can  be  approximately 
modeled  by  using  a  sequence  of  discrete  explosions,  separated  in  time  and  space,  followed 
by  superposition  of  the  single-source  surface  effects.  Further,  many  large  mining  explosions 
not  only  produce  large  surface  velocities,  but  also  cause  fracturing  and  eject  gas  and  dust 
into  the  atmosphere.  These  ejecta  produce  mass  coupling  to  the  atmosphere  and  subsequent 
air  flows  and  waves  as  well. 

Underground  nuclear  test  explosions,  which  are  usually  deeply  buried  so  as  to  avoid 
surface  fracturing  and  gas  venting,  produce  a  nearly  spherical  shock  wave  with  a  surface 
interaction  that  is  roughly  symmetric.  The  compressional  wave  from  this  source  produces  an 
upward  vertical  motion  followed  by  a  downward  motion.  Depending  on  the  size  of  the 
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nuclear  yield,  its  depth  and  the  degree  of  coupling  with  the  surrounding  rock,  the  surface 
boundary  can  experience  accelerations  well  over  lg.  Because  of  the  amplitude  of  the 
acceleration  at  the  surface,  the  surface  movement  is  usually  complicated  by  the  occurrence 
of  spallation  (separation  of  the  rock  along  bedding  planes),  so  that  several  peaks  in  the 
ground  acceleration  occur.  In  any  case,  the  peak  vertical  velocity  of  the  surface  directly 
above  the  source  can  vary  from  hundreds  of  meters  per  second,  for  the  largest  tamped 
explosions,  down  to  less  than  a  meter  per  second,  for  small  decoupled  explosions  detonated 
at  moderate  depths. 

Since  our  purpose  is  to  investigate  the  major  consequences  of  the  boundary  movements 
produced  by  these  sources,  we  will  focus  on  estimating  the  amplitudes  and  time-dependent 
wave-form  characteristics  of  the  gravity  waves  that  may  be  excited.  In  particular,  we  seek  to 
determine  whether  the  predicted  amplitudes  and  wave  characteristics  imply  effects  that  are 
easily  measurable  and  can  be  used,  along  with  seismic  signal  data,  to  distinguish  between 
different  source  types,  particularly  between  the  different  explosion  types. 

Qualitative  differences  in  the  relative  amplitudes  of  seismic  and  gravity  waves  produced 
by  the  different  seismic  sources  are  indicated  in  Table  1.  The  seismic  wave  amplitude 
differences  listed  are  based  on  theoretical  and  observational  results  (e.g.  EVERNDEN  etal., 
1987),  while  those  for  the  gravity  wave  amplitudes  are  based  on  inferences  from  the  limited 
observations  and  deductions  from  first  order  theoretical  considerations.  For  example,  in 
Table  1  we  infer  small  to  moderate  amplitude  gravity  wave  signals,  in  the  5  to  6  minute 
period  range,  for  thrust  and  normal  earthquakes.  This  inference  is  based  on  the  GPS 
measurements,  obtained  by  Calis  and  Minster  (1995),  following  the  rather  large  Northridge 
Earthquake.  On  the  other  hand,  only  very  small  gravity  wave  signals  are  expected  for  strike- 
slip  earthquakes  because  of  the  small  vertical  movements  to  be  expected  at  the  atmosphere 
boundary.  Likewise,  as  noted  earlier,  large  atmospheric  effects  have  been  observed  for  near 
surface  chemical  explosions  while  the  atmospheric  effects  from  more  deeply  buried 
explosions,  such  as  underground  nuclear  tests,  will  be  reduced  through  propagation  of  signals 
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to  the  surface,  as  well  as  by  the  possible  partial  decoupling  of  the  explosion  from  the 
medium  around  the  source  origin,  including  that  achieved  through  use  of  a  large  cavity  at 
the  detonation  point.  More  quantitative  and  precise  results  than  those  given  in  Table  1 
should  be  obtained  from  the  successful  application  of  a  non-linear  numerical  modeling 
approach  and  a  well  defined  observational  program.  Therefore,  at  this  stage,  we  hope  to 
establish  a  more  quantitative  basis  for  source  identification  through  a  well-developed 
predictive  capability. 

In  regard  to  modeling  the  atmospheric  excitation  by  seismic  sources,  we  only  require 
the  specification  of  the  particle  velocity  produced  by  seismic  waves  at  the  earth-atmosphere 
interface  in  order  to  express  the  coupling  of  the  seismic  source  to  atmospheric  flow  and  wave 
excitation  via  eq:  (15).  The  general  form  of  the  seismic  particle  velocity  at  the  earth's 
surface  due  to  a  shallow  seismic  source  is,  to  first  order,  described  by  a  maximum  velocity 
level  near  the  source  epicenter,  and  along  the  strike  of  the  fault  in  the  case  of  an  earthquake, 
with  rapid,  near  exponential,  fall-off  in  directions  away  from  the  fault  line  and  the 
epicenteral  zone.  The  time  history  of  the  normal  component  of  the  boundary  velocity  can 
also  be  approximated  by  a  simple  functional  variation,  where  we  expect  either  upward 
movement  followed  by  a  downward  rebound,  or  the  reverse  in  the  case  of  some  earthquake 
sources.  Of  course,  in  all  cases  there  will  be  local  deviations  from  this  uniform,  first  order, 
boundary  velocity  variation  due  to  medium  heterogeneities  and  anisotropic  characteristics, 
as  well  as  from  surface  asymmetries  and  non-linearities.  However,  these  variations  are  of 
short-wavelength  spatially  and  of  high  frequency  temporally.  They  are  therefore  averaged 
out  during  the  excitation  of  the  relatively  low  frequency,  long  wavelength  gravity  waves;  so 
that  in  modeling  gravity  wave  excitation  the  first  order  long  wavelength  approximation  of  the 
lower  atmosphere  boundary  velocity  movement  is  adequate. 

The  previous  considerations  suggest  that  a  good  first  order  approximation  to  the 
boundary  velocity  can  be  obtained  by  taking  the  normal  component  of  the  seismic  particle 
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velocity  to  be  of  Gaussian  form  on  the  surface  (S0)  which  we  take  to  be  the  plane  X3  =  0. 
In  particular,  on  S0: 
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withvn  denoting  the  solid  medium  particle  velocity  normal  to  S0  and  x.,  i  =  1,2  and  t 


being  the  epicenteral  coordinates  of  the  source  and  the  origin  time,  respectively.  Here  also 
Ax,  (with  i  =  1,2)  and  At  are  the  Gaussian  scale  factors  for  the  spatial  and  temporal  time 
variations.  The  parameter  D0  is  the  (positive)  magnitude  of  the  maximum  dynamic  surface 
displacement.  Taking  the  time  derivative  and  applying  the  continuity  condition  for  the 
particle  velocities  across  S0,  with  w  denoting  the  vertical  component  of  the  particle  velocity 
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This  expression  contains  enough  flexibility  to  provide  first  order  representations  of  the 
boundary  velocity  in  the  near  field  for  all  the  seismic  sources  of  interest,  since  the  spatial 
scale  factors  Ax,  and  Ax2  can  be  chosen  differently  and  the  surface  velocity  variations  for 
the  different  source  geometries  can  all  be  well  represented. 

Typical  values  of  the  parameters  Ax,  for  small  to  moderate  earthquakes  are  in  the  range 
from  about  1  to  5  km.  These  choices  for  the  spatial  parameters  are  roughly  appropriate  to 
earthquakes  with  magnitudes  in  the  range  from  about  3  to  5.  The  time  parameter  At  is 
proportional  to  the  dominant  period  of  the  particle  velocity  variation  at  the  surface; 
specifically  it  is  about  one-fourth  of  the  period.  In  this  near  field  distance  range  an 
appropriate  range  for  At  can  be  obtained  from  observations  of  the  dominant  frequency  of 
near  field  seismic  waves  for  shallow  sources,  with  depths  less  than  5  km.  For  explosions  in 
the  3  to  5  magnitude  range  the  dominant  frequencies  are  from  about  25  to  5  HZ,  respectively 
(EVERENDEN  etal.,  1987).  Thus,  the  appropriate  At  range  is  from  about  0.1  to  .05.  For 
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earthquakes  in  the  same  magnitude  range  the  dominant  seismic  signal  frequencies  are 
somewhat  lower,  so  that  the  appropriate  range  for  At  is  from  about  .05  to  .25  sec. 

Based  on  observations  of  accelerations  greater  than  1  g  at  the  surface  for  shallow 
moderate  sized  earthquakes,  along  with  the  associated  measurements  of  peak  velocities  that 
range  up  to  several  tens  of  meters  per  second,  we  infer  that  the  displacement  parameter  D0 
appearing  in  (24)  should  have  a  range  of  from  about  .5  m  to  about  2.5  m  for  shallow 
earthquakes  with  magnitudes  ranging  from  3  to  5.  For  underground  explosions,  however,  the 
signal  frequencies  are  higher  and  D0  values  are  lower  (.1  m  to  .5  m)  for  this  magnitude 
range.  These  values  of  D  0  are  also  consistent  with  observations. 

Therefore,  the  first  order  representation  in  (24)  using  the  various  parametric  values 
indicated  above,  can  be  made  consistent  with  average  near-field  seismic  observations  from 
small  to  moderate  earthquakes  and  exlosions  and  should  provide  good  first  order  estimates 
for  the  boundary  velocity  variations.  The  examples  given  in  the  next  section  will  be 
representative  of  the  excitation  of  low  frequency  atmospheric  gravity  waves  by  seismic 
events  in  the  interesting  magnitude  range  from  3  to  5.  For  this  purpose  we  will  use  an 
earthquake  near  magnitude  4  and  use  a  maximum  value  of  about  4  m/sec  for  w,  with  At  =  .1 
sec.  (so  D0  =  .2  m)  and  with  Axx  =  2Ax2  =  4  km.  These  parameters  are  roughly  appropriate 
for  an  earthquake  of  magnitude  4  at  a  depth  of  about  5  km.  For  shallower  depths  the  values 
of  w  and  D0  increase.  At  a  depth  of  1  km  the  maximum  value  for  w  would  be  near  50  m/sec, 
with  D0  about  2  m. 

Modeling  of  Gravity  Waves  and  Associated  Ionospheric  Electron  Density  Variations 

The  boundary  source  represented  in  (24),  with  parameters  as  enumerated  for  a 
magnitude  4  earthquake,  produces  a  compressive  movement  in  the  atmosphere  followed  by 
an  expansion  when  the  event  has  a  thrust  mechanism.  The  initial  increased  density  from  the 
compressive  pulse  in  the  atmosphere  is  propagated  upward  more  slowly  than  the  following 
dilatational  pulse  which  has  increased  buoyancy  due  to  reduced  density.  One  also  finds  that 
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the  time-dependent  transient  pulse  in  3-dimensions  has  nonsymmetric  flows  that  are  critical 
to  the  upward  motion  of  the  buoyancy  wave.  In  particular,  as  a  function  of  time  the 
transients  propagate  upward  with  advective  motion  occurring  both  vertically  and  horizontally 
and  a  sequence  of  spiral  circulation  patterns  develops  through  asymmetric  screw  modes 
across  the  horizontal  cross-section,  as  previously  argued  by  NYGREN  et  al.,  (1984).  The 
circulation  patterns  across  the  vertical  cross-section  are  characterized  by  upward  central 
motions  of  the  lighter  matter,  which,  at  the  neutral  buoyancy  level,  pushes  outward  to  the 
side.  Then,  as  the  pulse  centroid  moves  through  this  level  the  air  mass  is  pulled  back  toward 
the  center.  The  advected  air  mass  therefore  tries  to  remain  in  its  horizontal  stratification  in 
order  to  minimize  changes  in  its  gravitational  potential.  Further,  energy  and  momenta  are 
circulated  through  the  plume  by  traveling  waves  through  the  circulation  pattern,  as  was  also 
observed  by  WEIL  (1988).  These  effects  have  been  observed  in  the  real  atmosphere  when 
thermals  and  plumes  propagate  upward  with  similar  asymmetric  shearing  flows  (KUETTNER 
etal.,  1987).  We  also  observe  that  the  pulse  increases  its  horizontal  wavelengths  as  time 
progresses,  with  800  km  wide  grids  insufficient  to  map  the  whole  pulse  structure  at  an 
altitude  of  100  km. 

Figure  2  shows  the  normalized  horizontal  velocity  as  a  function  of  time  and  location 
along  an  axis  through  the  grid,  for  this  seismic  source.  The  altitude  of  the  flow  is  about  170 
km  and  amplitude  values  are  normalized  to  the  velocity  of  sound.  The  flipping  of  the 
velocity  direction  in  time  is  regular,  with  the  change  in  direction  along  the  axis  occurring 
directly  over  the  source.  The  flipping  of  the  direction  of  flow  velocity  at  sharp  density 
gradients  was  predicted  by  the  linear  theory  approximations  for  gravity  waves  (HINES, 
1960,1974),  and  was  confirmed  numerically  for  two  dimensional  systems  by  GREENE  and 
WHITAKER  (1968)  who  found  gravity  waves  forming  at  the  120  km  altitude.  Thus,  it 
appears  as  if  a  source  of  thermospheric  modes  exists  at  this  altitude  (FRANCIS,  1975).  The 
downward  slanted  direction  of  the  phase  propagation  is  readily  observed,  as  predicted  by  the 
linear  approach  of  HINES  (1960).  Figure  3  shows  the  behavior  of  the  normalized  perturbed 
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Figure  2.  Normalized  horizontal  velocity,  scaled  by  the  surface  air  velocity  of  3.4  m/sec,  for 
a  gravity  wave  as  a  function  of  time  after  its  initial  onset.  The'  x  component  of  the  particle 
velocity  at  a  height  of  170  km  above  the  surface  is  shown.  The  boundary  source  has  a 
Gaussian  spatial  distribution,  with  Ax,  =2Ax2  =  4  km.  The  boundary  velocity  time  variation 

corresponds  to  the  time  derivative  of  Gaussian  variation,  with  At  =  .1  sec.  The  maximum 
boundary  velocity  value  is  4  m/sec.  A  dynamic  drag  coefficient  of  .01  for  the  atmosphere- 

ionosphere  is  assumed  along  with  a  choice  of  a  numerical  time  step  of  .1  sec  in  the  finite 
difference  computation. 
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Figure  3.  Normalized  pressure  (scaled  by  atmospheric  pressure  at  the  earth's  surface)  of  a 
gravity  wave  as  a  function  of  time  after  initial  onset.  The  source  type  and  drag  coefficient 
are  identical  to  those  for  Figure  2,  and  the  height  is  also  at  170  km. 
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pressure  field  under  the  same  conditions  as  in  Figure  2,  with  all  these  results  being 
consistent  with  the  earlier  predictions  from  2  D  numerical  modeling  and  the  linear  analytical 
estimates. 

Figure  4  shows  the  temperature  disturbance  above  the  source  as  a  function  of  time  at 
different  heights.  An  initial  large  pulse  decays  in  amplitude  with  oscillations  at  about  a 
5  minute  period.  The  gravity  wave  interacting  with  the  steep  gradients  at  the  thermopause  is 
thought  to  cause  these  leaky  modes  and  their  periods  are  seen  to  vary  from  3  to  5  minutes, 
and  includes  the  lowest  orders  of  acoustic  and  gravity  modes.  Such  oscillations  have  been 
recorded  in  numerous  observations  of  waves  produced  by  near  surface  sources  (GEORGES, 
1968)  and  in  observations  of  ionosondes  near  Space  Shuttle  ascents  (JACOBSON  and 
CARLOS,  1994). 

To  first  order,  the  electrons  in  the  ionosphere  are  assumed  to  move  with  the  flow  of  the 
dominant  neutrals  and  the  change  in  the  electron  density  is  calculated  from  the  conservation 
law  in  equ:  (17),  where  integration  over  time  gives  the  total  electron  density  variation.  For 
the  boundary  velocity  level  used  for  the  magnitude  4  earthquake  we  find  changes  in  electron 
density  that  are  somewhat  less  than  one  tenth  of  a  percent.  When  the  drag  and  dissipative 
effects  are  small,  we  obtain  somewhat  larger  peak  amplitudes  and  also  find  oscillations  in 
the  electron  density  changes  that  are  induced  by  flipping  of  the  gravity  wave  at  the 
thermocline. 

Figure  5  shows  the  electron  density  changes  in  time  at  various  altitudes  between  125 
km  and  250  km  for  the  representative  magnitude  4  earthquake  source.  These  results  show 
the  characteristic  oscillations  imposed  on  the  temporally  and  spatially  complex  pulse  when 
the  turbulent  drag  coefficient  is  small.  When  the  dissipation  is  increased,  these  oscillations 
are  not  apparent.  Further,  the  dissipative  mechanisms  can  be  adjusted  to  vary  the  decay  of 
the  pulse  in  time  and  with  height. 

Figure  6  compares  the  electron  density  variations  in  the  ionosphere  due  to  gravity  waves 
from  a  source  of  the  form  used  previously,  with  the  same  pulse  duration  but  with  different 
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DEPTH  below  250km  (5km  units) 


Figure  5.  Electron  density  changes  (normalized  by  1012/m3)  after  the  onset  of  a  gravity 
wave,  for  altitudes  between  125  km  and  250  km.  The  boundary  velocity  source  is  as  in 
Figure  2.  The  dynamic  drag  coefficient  is  .001  in  this  example. 


ELECTRON  DENSITY  CHANGE  (normalised) 


Figure  6.  Normalized  electron-density  changes  due  to  a  gravity  wave  as  a  function  of  time 
after  initial  onset.  Depth  below  250  km  is  indicated  on  the  third  axis.  The  boundary  velocity 
source  has  the  same  form  as  described  in  Figure  2,  except  that  maximum  velocity  values  are 
varied.  The  drag  coefficient  was  taken  to  be  .02.  Case  (A):  The  scale  velocity  for  the 
source  is  3  m/sec.  Case  (B):  The  source  has  a  scale  velocity  of  15  m/sec.  (The  actual 
electron  density  change  is  10  times  the  value  on  the  scale  for  case  B). 
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maximum  boundary  velocities  appropriate  to  a  shallower  event.  These  simulations  assumed 
a  high  dynamic  drag  coefficient  of  .02.  The  results  have  qualitative  similarity  in  magnitude 
to  observed  electron  density  changes  above  near-surface  sources  (PIERCE  et  al.,  1971, 
GOKHBERG,  1991)  but  lack  oscillatory  character  due  to  the  high  drag  coefficient  assumed. 
Of  course,  if  the  drag  coefficient  is  reduced  oscillatory  variations  like  that  in  Figure  5  are 
obtained,  with  amplitudes  somewhat  above  those  shown  in  Figure  6. 

Comparisons  of  Modeling  Results  with  Observed  Electron  Density  Variations 
Inferred  from  GPS  Data 

On  January  17,  1994,  at  4:31  am  local  time,  the  Northridge  earthquake  occurred  in 
southern  California  about  30  km  north  of  the  Los  Angeles  basin  (HAUKSSON  et  al.,  1994). 
This  magnitude  6.6  earthquake  ruptured  with  a  thrust  mechanism  of  a  10  km  long  fault 
segment  striking  N120E  and  dipping  45°  to  the  southwest  (THIO  and  KANAMORI,  1994). 
The  rupture  started  at  a  19  km  depth  and  propagated  upward  to  a  depth  of  about  5  km,  but 
did  not  reach  the  surface  (WALD  and  HEATON,  1994).  This  produced  40  cm  of  static 
vertical  displacement  directly  above  the  fault  plane  (HUDNUT  et  al,  1994).  Ground 
accelerations  of  more  than  1.5  g  were  recorded  5  km  from  the  epicenter  (EGAL  etal.,  1994). 

The  Northridge  epicenter  was  located  within  the  permanent  network  of  GPS  stations 
operating  in  southern  California  (PGGA  network,  BOCK  et  al.,  1993)  and  CALAIS  and 
MINSTER  (1995)  used  the  GPS  data  from  the  PGGA  network  to  estimate  fluctuations  of  the 
ionospheric  electron  content  at  distances  up  to  1000  km  from  the  source.  Specifically,  since 
the  radio  signals  broadcast  at  1.57542  GHz  (fi)  and  1.2276  GHz  (f2)  by  GPS  satellites  are 
dispersively  delayed  along  their  path  by  interactions  with  free  electrons  in  the  ionosphere,  so 
that  the  differential  delay  between  the  fi  and  f2  frequencies  is  proportional  to  the  electron 
density  along  the  ray  path  (KLOCUBHAR,  1985),  then  measured  phase  delays  can  be  used 
to  calculate  time  dependent  electron  density  variations.  Therefore,  using  phase  and  P-code 
delays  recorded  on  the  two  GPS  frequencies,  CALAIS  and  MINSTER  made  direct 
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measurements  of  the  Total  Electron  Content  (TEC)  of  the  ionosphere,  corresponding  to  the 
electron  density  integrated  along  the  line  of  sight. 

The  TEC  values  along  the  signal  paths  at  the  GPS  receiving  stations  were  computed 
and  then  corrected  to  a  "vertical  electron  content"  (VEC)  above  the  receivers  in  order  to 
remove  the  geometric  effect  of  the  different  satellite  paths  to  the  receivers.  Since  the  power 
spectrum  of  the  VEC  variation  is  dominated  by  very  low  frequency  components  introduced 
by  the  solar  cycle,  it  is  necessary  to  filter  out  this  large  noise  component  in  order  to  extract  a 
signal  corresponding  to  the  electron  density  variations  introduced  by  any  gravity  wave. 
Therefore,  CALAIS  and  MINSTER  band-pass  filtered  the  VEC  data  in  the  period  range  from 
about  3  to  10  minutes,  based  on  the  modeling  predictions  shown  in  Figure  5. 

The  results  of  applying  the  simple  band  pass  filtering  to  the  GPS  data  at  several 
receiver-satellite  pairs  near  the  earthquake  epicenter  (within  200  km  or  so)  are  shown  in 
Figures  7  and  8.  In  Figure  7  a  single  (typical)  recording  of  the  VEC  by  a  ground  GPS 
receiver,  covering  the  time  frame  from  just  before  the  earthquake  until  just  after,  is  shown. 
A  signal  arriving  shortly  after  the  earthquake,  with  a  signal  to  noise  ratio  of  about  2  or 
slightly  greater,  can  be  seen  in  this  record.  This  signal  has  a  period  of  about  5  minutes  and  a 
time  variation  approximating  that  predicted  in  Figure  5.  Naturally,  because  of  the  movement 
of  the  satellite  there  is  some  distortion  of  the  VEC  signal  observed.  Nevertheless,  the 
agreement  between  prediction  and  observation  is  quite  good  in  terms  of  wave  form  similarity 
and  the  predicted  and  observed  dominant  period. 

Since  there  are  many  sources  of  ionospheric  gravity  waves,  it  could  only  be  fortuitous 
that  a  larger  signal  above  the  noise  level  is  seen  shortly  after  the  earthquake  at  a  single  GPS 
receiver.  However,  Figure  8  shows  several  receiver  recordings,  band  pass  filtered  as  in 
Figure  7,  that  all  show  a  signal  above  the  background  level  arriving  at  different  times  after 
the  earthquake.  In  fact,  since  the  records  are  aixanged  with  distance  of  the  receiver-satellite 
signal  path  from  the  event  epicenter,  which  increases  from  the  bottom  to  the  top  of  the 
record  section  shown,  it  can  be  seen  that  the  VEC  signal  has  a  move  out  (variable  time 
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arbitrary  units  (same  for  aU  traces) 


zero-phase  filtering 


Figure  7  Pass-Band  Filtered  (3-10  minute  period  range)  VEC  time  recording  for  a  single 
satellite-receiver  pair  (site  MATH,  satellite  13).  The  dashed  line  indicates  the  ori¬ 
gin  time  of  the  earthquake  (magnitude  6.6,  ).  A  signal,  above  background 

by  about  a  factor  of  2,  is  detected  with  onset  between  4:30  to  5:00  local  time. 
(Figure  from  Calais  and  Minster,  1994) 


pass-band  filtered  10-3  inn 


Figure  8  Observed  filtered  VEC  recordings  at  a  number  of  receiver-satellite  pairs.  The 
recordings  are  ordered  by  increasing  distance  of  the  signal  path  from  the  epi¬ 
center,  with  the  closest  at  the  bottom.  The  dotted  line  shows  the  earthquake  ori¬ 
gin  time  and  the  numbers  to  the  right  indicate  receiver-satellite  pairs  for  each 
recording.  All  the  recordings  show  evidence  of  a  signal,  with  a  time  delay  of  its 
onset  increasing  with  increasing  distance  of  the  receiver-satellite  path  from  the 
event  epicenter,  as  would-be  expected  for  a  slowly  moving  ionospheric  gravity 
wave  perturbation  of  the  electron  density.  (Figure  from  Calais  and  Minster,  1994) 
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delay)  that  increases  with  distance.  This,  of  course,  would  be  expected  for  an  ionospheric 
electron  density  perturbation  arising  from  a  gravity  wave  source  near  the  epicenter  of  the 
earthquake.  The  propagation  speed  indicated  by  this  sequence  of  signal  time  delays  is 
between  300  and  600  m/sec,  a  range  consistent  with  gravity  wave  speeds  in  the  ionosphere. 
Furthermore,  the  initiation  of  the  disturbance  has  a  delay  time  consistent  with  the  earthquake 
excitation  of  gravity  waves  near  the  epicenter  at  the  time  of  the  event.  Therefore,  there  is 
little  doubt  that  the  perturbations  in  the  electron  density  along  all  the  signal  paths  are  due  to 
propagating  gravity  waves  which  excited  by  this  earthquake. 

A  somewhat  more  complete  comparison  of  the  observations  and  the  modeling 
predictions  can  be  made  using  previous  modeling  results.  Of  course  the  magnitude 
4  earthquake  assumed  in  the  computations  for  the  results  in  Figure  5  is  quite  different  than 
the  Northridge  earthquake,  which  was  of  larger  magnitude  (6.6)  and  was  deeper,  with  a 
hypocentral  depth  of  about  19  km.  However,  adjusting  for  the  depth  differences  as  well  as 
the  effects  of  the  magnitude  differences  between  the  events,  indicates  that  the  Northridge 
event  would  have  a  boundary  velocity  about  as  large  as  that  for  the  shallower  "standard" 
earthquake,  but  with  values  of  Ax;  somewhat  larger  than  those  used  for  the  representative 
magnitude  4  earthquake.  Thus,  the  Northridge  earthquake  would  be  expected  to  produce 
somewhat  larger  atmospheric  gravity  waves,  but  of  the  same  order  of  magnitude.  This 
estimate,  in  fact,  is  consistent  with  the  maximum  VEC  predicted  by  the  results  in  Figure  5. 
That  is,  the  maximum  unnormalized  electron  density  variation  from  the  predictions  is 
approximately  (3  x  10"4  )  x  1012  el/m3,  or  3  x  108  el/m3.  Converted  to  a  VEC  variation  by 
simply  multiplying  by  the  thickness  (h)  of  the  ionosphere,  gives  about  3  x  10 14  el/m2.  The 
value  observed  for  the  Northridge  earthquake  is  (coincidentally)  very  close  to  this  value. 

Of  course  these  comparisons  are  quite  crude,  since  there  has  been  no  systematic  effort 
to  vary  atmospheric  dissipation  in  the  modeling  in  order  to  match  the  dominant  period  of  the 
observed  gravity  waves,  and  the  source  differences  are  only  roughly  estimated. 
Nevertheless,  the  rough  agreement  with  the  predictions  from  a  "standard"  earthquake  is 
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significant,  and  there  can  be  little  doubt  that  relatively  minor  adjustments  to  the 
atmospheric-ionospheric  dissipation  parameters,  as  well  as  to  source  boundary  velocity 
parameters,  could  produce  quite  detailed  fits  to  the  observations. 

SUMMARY  AND  CONCLUSIONS 

Non-linear  numerical  modeling  in  3-D  has  been  used  to  investigate  the  excitation  of 
atmospheric  and  ionospheric  disturbances  due  to  shallow  (upper  crustal)  seismic  sources. 
Grid  zones  as  large  as  800  by  800  kilometers  in  area  and  400  km  n  altitude  were  used  with 
grid  dimensions  between  2  and  15  km  also  used.  The  modeling  employs  non-local 
differencing  for  the  non-linear  terms  in  the  equations  and  results  in  a  non-linear  integration 
algorithm  that  is  found  to  be  essential  for  accuracy  and  stability.  Altitude  dependent  wave 
dissipation  effects,  such  as  arise  from  interaction  of  the  gravity  waves  with  background  wind 
turbulence,  are  modeled  using  a  drag  force  proportional  to  the  negative  of  the  wave  particle 
velocity.  Predictions  of  gravity  wave  characteristics  are  found  to  be  very  sensitive  to  the 
magnitude  of  the  dynamic  drag  coefficient,  with  oscillatory  gravity  waves  occurring  all  the 
way  into  the  ionosphere  for  values  below  .01.  For  a  larger  average  drag  coefficient  the 
waves  are  over  damped  and  oscillations  are  absent. 

Results  from  numerical  simulations  agree  with  predictions  from  linear  theory 
approximations  and  also  show  non-linear  characteristics  that  have  been  observed 
experimentally.  A  comparison  of  model  predictions  of  electron  density  fluctuations  due  to 
earthquake  generated  atmospheric  gravity  waves  with  observations  of  the  variations  in  the 
Vertical  Electron  Concentration  (VEC)  obtained  from  GPS  receivers  after  the  Northridge 
earthquake  in  California  showed  agreement  in  the  oscillatory  wave  form  and  dominant  (near 
5  minutes)  period  of  the  oscillations.  In  addition,  the  amplitude  levels  of  the  observations 
and  predictions  were  consistent  with  each  other.  The  agreement  with  observations  implies 
that  a  low  value  for  the  average  dynamic  drag  coefficient  is  appropriate  and  that  small  to 
moderate  shallow  earthquakes,  of  thrust  or  normal  type,  will  produce  gravity  waves  of 
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sufficient  amplitude  to  be  detectable  through  the  perturbations  they  cause  in  the  ionospheric 
electron  density.  Additionally,  GPS  satellites  and  receivers  provide  a  ready  made  network 
for  such  measurements  using  the  approach  followed  by  CALAIS  and  MINSTER  (1995)  in 
their  study  of  the  signals  from  the  Northridge  earthquake. 

The  opportunity  to  use  GPS  observations  and  the  non-linear  modeling  capability  to  study 
atmosphere  and  ionosphere  properties  seems  evident  since  there  would  appear  to  be  an 
abundance  o  shallow  seismic  sources  which  can  excite  gravity  waves  producing  measurable 
effects.  In  particular  shallow,  upper  crustal  thrust  and  normal  earthquakes,  volcanic 
eruptions,  large  subsurface  (chemical  and  nuclear)  explosions,  as  well  as  much  smaller 
chemical  explosions  at  the  surface,  should  all  produce  measurable  effects.  Further,  it  should 
be  possible  to  use  the  electromagnetic  sounding  measurements,  along  with  seismic 
measurements,  to  differentiate  between  these  various  seismic  sources.  Therefore,  a  dual 
field  sensing  approach,  involving  GPS  and  seismic  field  sensors,  could  be  used  to  better 
distinguish  between  large  industrial  explosions  and  underground  nuclear  tests,  particularly 
tests  producing  signals  that  have  been  reduced  by  decoupling  techniques. 
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